
%% load in previously selected referencing information:
[fn,fd] = uigetfile('*.mat','select referencing information .mat file');
cd(fd);
load(fn);

%% load in vessel mask:
load(vesselmaskname,'vesselMasks','avios','refim_windowmask')     %avios is what we're going to be registering everything to; there is ONE of these per mouse 
refim = avios;              %need to make a reference ios image for each session

%% load in background:
cd(bgdirectory);
bg_ios = loadRawStack_v1(rootname,localimagesize);

load(bgmetaname,'movindicator','rectPositions','framesPerStim')

bg_ios = reshape(bg_ios,size(bg_ios,1),size(bg_ios,2),framesPerStim,length(movindicator));

bguid = [reshape(rectPositions,2,size(rectPositions,3));movindicator];
[uid_bg,~,bgID_idx] = unique(bguid','rows');


%% load in metadata

load(metaname,'movindicator','rectPositions','framesPerStim','preStimTime','totalFrameTime','stimparams');
stimuid = [reshape(rectPositions,2,size(rectPositions,3));movindicator];
[uid,~,stimID_idx] = unique(stimuid','rows');
totalFrameTime = totalFrameTime/1000;       %convert into seconds

time = (1:(framesPerStim))*totalFrameTime;
prestimframes = 1:sum(time<preStimTime);
stimframes = (prestimframes(end)+1):sum(time<(preStimTime + stimparams{1}.totalStimLength));
iosResponseFrames = find(time>(preStimTime+1),1):find(time>(preStimTime + stimparams{1}.totalStimLength+1),1);
%% will want to get ONE very good tranformation matrix to map the data from the current session to the master session
% to do this you need to load in some ios data from the current session,
% make an average image; do NOT need to do this if this is the reference
% session:
if ~isReferenceSession
    cd(movdirectory);
    curidx = 1:(max(prestimframes));
    localavios = loadRawStack_v1(rootname,localimagesize,curidx);
    localavios = mean(localavios,3);
    A = localavios;A(~windowmask) = median(A(:));
    B = refim; B(~refim_windowmask) = median(B(:));
    sessionRegister = registerMultipleWidefieldImages_v1_2(A,B);
    localavios_transformed = imwarp(localavios,sessionRegister.Transformation,'OutputView',imref2d([size(refim,1),size(refim,2)]));
    figure; imshowpair(refim,localavios_transformed,"checkerboard");
    clear A B
else
    localavios = avios;
end
%% loop through movies and do the tracking (for .raw formatted image sequences):
cd(movdirectory);
totalframes = framesPerStim*length(movindicator);
traj = [];
gcampmov_full = [];
parenchymaMask = false(size(bg_ios,1),size(bg_ios,2),length(movindicator));
roughvesselmask = false(size(bg_ios,1),size(bg_ios,2),length(movindicator));

avgParenchymaIOS_normalized = zeros(size(parenchymaMask));
for n = 1:length(movindicator)
    curidx = ((n-1)*framesPerStim + 1):(n*framesPerStim);
    iosmov = loadRawStack_v1(rootname,localimagesize,curidx);
    curbg = bgID_idx==stimID_idx(n);
    iosmov = iosmov-bg_ios(:,:,:,curbg);
    if ~isReferenceSession
        [iosmov,curtraj] = trackOneWidefieldMovie_iosonly(iosmov,windowmask,localavios,refim,vesselMasks, sessionRegister);
    else
        [iosmov,curtraj] = trackOneWidefieldMovie_iosonly(iosmov,windowmask,localavios,[],vesselMasks);
    end
    [parenchymaMask(:,:,n),roughvesselmask(:,:,n)] = getParenchymaPixels(mean(iosmov,3),windowmask);
    
    for n2 = 1:size(iosmov,3)
    
        temp = iosmov(:,:,n2);
        temp(~parenchymaMask(:,:,n)) = nan;
        temp(~windowmask) = 0;
        iosmov(:,:,n2) = temp;
    end
    temp = mean(iosmov(:,:,prestimframes),3,'omitnan');
    temp = iosmov./repmat(temp,1,1,size(iosmov,3));
    avgParenchymaIOS_normalized(:,:,n) = mean(temp(:,:,iosResponseFrames),3,'omitnan');

    traj = [traj,curtraj];
    disp(['tracked movie # ',num2str(n)])
end
numTrackedSegs = size(traj(1).vesselCent,1);
disp('FINISHED TRACKING DILATIONS')
%% calculate baseline levels for vessel total intensity:
clear iosmov mov curgcampmov curtraj
baselineVesselInten = zeros(numTrackedSegs,1);
for n = 1:length(traj)
    baselineVesselInten = baselineVesselInten+mean(traj(n).vesselTotalInten(:,prestimframes),2,'omitnan');
end
baselineVesselInten = baselineVesselInten/length(traj);

%% calculate  baseline levels for surround total intensity:
baselineSurroundInten = zeros(numTrackedSegs,1);
for n = 1:length(traj)
    baselineSurroundInten = baselineSurroundInten+mean(traj(n).surroundInten(:,prestimframes),2,'omitnan');
end
baselineSurroundInten = baselineSurroundInten/length(traj);

    
%% run processing on gcamp movie AFTER the fact to allow for global baseline correction:
clear iosmov curgcampmov curtraj

%% save data so it's not lost
cd(savedirectory)
save(savename,'baselineVesselInten','traj','stimID_idx','uid','parenchymaMask','roughvesselmask','baselineSurroundInten','avgParenchymaIOS_normalized');


